rm(list=ls(all=TRUE))

library(foreign)
library(MASS)
library(memisc)
library(plyr)
library(stringr)
library(utils)
library(stats)
library(pscl)
library(MuMIn)
library(BMA)
library(car)

if(Sys.info()["user"]=="IV")
  setwd("/Users/IV/")
if(Sys.info()["user"]=="benjaminbarber")
  setwd("/Users/IV/")


Combined_Data_Pop <- read.dta("/Users/IV/Dropbox/Nazi/Data/Master Data for Analysis.dta")





####### Setting Up Variable Lists ######

must_haves <- "as.factor(muster_year)"
region_fixed_effects <- "in_welfare_per1000 + war_per1000 + sozialrentner_per1000 + bigd + logtaxprop + c25juden_share + c25kath_share + c33erlos_share + c33brlos_share + c25anges_share + c25arbei_share"
population <- "population_fix + I(population_fix^2) + I(population_fix^3) + I(population_fix^4) + I(population_fix^5)"
sig_str <- "avg_max_signal_strength_fix"
logit_sig_str <- "avg_max_signal_strength_fix_log"


choice <- c("catholic", "wounded", "married", "class", "hum_cap", "Relative_Height", "Relative_Weight", "num_siblings", "num_children", "nwounds", "NSDAP_Fix", "nazi2", "urban", "age_at_muster_relative18")
choice.n <- length(choice)


num_of_models <- 0
for(i in 1:choice.n){
	X <- factorial(choice.n)/(factorial(i)*factorial(choice.n - i))
	num_of_models <- num_of_models + X
}
num_of_models

Combined_Data_Pop$num_children <- ifelse(is.na(Combined_Data_Pop$num_children)==TRUE,0, Combined_Data_Pop$num_children)

N <- 0

Sala_i_Martin_unweight <- matrix(NA, ncol=4, nrow=3)
Sala_i_Martin_weight <- matrix(NA, ncol=4, nrow=3)


############# EBA Analysis #############


###### Medals ######

## DV: Decorations: decorated
## IV: Radio Signal Strength
N <- N+1

id <- unlist(
	lapply(1:choice.n,
		function(i)combn(1:choice.n,i,simplify=FALSE)
		)
	,recursive=FALSE)

formulas <- sapply(id, function(i)
	paste("dec ~ 1", sig_str, must_haves, population, region_fixed_effects, paste(choice[i],collapse=" + "), sep=" + ")
	)

#testing 
# m <- lm(as.formula(formulas[1]), data= Combined_Data_Pop)
	
savefunc <- function(i){
	m <- lm(as.formula(i), data= Combined_Data_Pop)
	f <- summary(m)$r.squared
	c(summary(m)$coef[2,], f)	
}

X <- lapply(formulas,savefunc)

point <- NULL
se <- NULL
z <- NULL
r2 <- NULL

for(i in 1:length(X)){
	point <- c(point, X[[i]][[1]])
	se <- c(se, X[[i]][[2]])	
	z <- c(z, X[[i]][[3]])
	r2 <- c(r2, X[[i]][[5]])		
}

lower <- point + qnorm(0.025)*se
upper <- point + qnorm(0.975)*se

dv <- data.frame(point, lower, upper, se, z, r2)
dv <- dv[order(dv$point),]
colnames(dv) <- c("point", "lower", "upper", "se", "z", "r2")
dv$number <- seq(1:dim(dv)[1])
dv$weight <- dv$r2/sum(dv$r2)

# Levine-Renelt
dv$lower[1]
dv$upper[dim(dv)[1]]

c(dv$point[1], dv$se[1], dv$z[1])
c(dv$point[dim(dv)[1]], dv$se[dim(dv)[1]], dv$z[dim(dv)[1]])

# Sala--i-Martin
# unweighted:
c(mean(dv$point), mean(dv$se), mean(dv$z), max(pnorm(0,mean(dv$point), mean(dv$se), lower.tail=T), pnorm(0,mean(dv$point), mean(dv$se), lower.tail=F)))

Sala_i_Martin_unweight[N,] <- c(mean(dv$point), mean(dv$se), mean(dv$z), max(pnorm(0,mean(dv$point), mean(dv$se), lower.tail=T), pnorm(0,mean(dv$point), mean(dv$se), lower.tail=F)))

# weighted
weight_point <- sum(dv$weight*dv$point)
weight_se <- sum(dv$weight*dv$se)
weight_z <- sum(dv$weight*dv$z)

c(weight_point, weight_se, weight_z, max(pnorm(0,weight_point, weight_se, lower.tail=T), pnorm(0,weight_point, weight_se, lower.tail=F)))
Sala_i_Martin_weight[N,] <- c(weight_point, weight_se, weight_z, max(pnorm(0,weight_point, weight_se, lower.tail=T), pnorm(0,weight_point, weight_se, lower.tail=F)))




# Graph
quartz(type="pdf", width=5, height=5, file="/Users/IV/Dropbox/Nazi/Output/EBA_Decorated_Signal_Strength_ALL.pdf")
plot(NULL, xlim=c(1, length(dv$number)), ylim=c(0, max(dv$upper)), ylab=expression(paste(beta, " Radio Signal Strength")), xlab="", main="EBA - Decorations")
for(i in 1:length(dv$number)){
	lines(c(i,i),c(dv$lower[i], dv$upper[i]), lwd=1.5, col="grey")
}
points(dv$number, dv$point, pch=16,cex=0.5)
abline(h=0, lwd=2, col="red")
dev.off()





## DV: High Decoration
## IV: Radio Signal Strength
N <- N+1

id <- unlist(
	lapply(1:choice.n,
		function(i)combn(1:choice.n,i,simplify=FALSE)
		)
	,recursive=FALSE)

formulas <- sapply(id, function(i)
	paste("High_Medal ~ 1", sig_str, must_haves, population, region_fixed_effects, paste(choice[i],collapse=" + "), sep=" + ")
	)
	
savefunc <- function(i){
	m <- lm(as.formula(i), data= Combined_Data_Pop)
	f <- summary(m)$r.squared
	c(summary(m)$coef[2,], f)	
}

X <- lapply(formulas,savefunc)

point <- NULL
se <- NULL
z <- NULL
r2 <- NULL

for(i in 1:length(X)){
	point <- c(point, X[[i]][[1]])
	se <- c(se, X[[i]][[2]])	
	z <- c(z, X[[i]][[3]])
	r2 <- c(r2, X[[i]][[5]])		
}

lower <- point + qnorm(0.025)*se
upper <- point + qnorm(0.975)*se

dv <- data.frame(point, lower, upper, se, z, r2)
dv <- dv[order(dv$point),]
colnames(dv) <- c("point", "lower", "upper", "se", "z", "r2")
dv$number <- seq(1:dim(dv)[1])
dv$weight <- dv$r2/sum(dv$r2)

# Levine-Renelt
dv$lower[1]
dv$upper[dim(dv)[1]]

c(dv$point[1], dv$se[1], dv$z[1])
c(dv$point[dim(dv)[1]], dv$se[dim(dv)[1]], dv$z[dim(dv)[1]])

# Sala--i-Martin
# unweighted:
c(mean(dv$point), mean(dv$se), mean(dv$z), max(pnorm(0,mean(dv$point), mean(dv$se), lower.tail=T), pnorm(0,mean(dv$point), mean(dv$se), lower.tail=F)))

Sala_i_Martin_unweight[N,] <- c(mean(dv$point), mean(dv$se), mean(dv$z), max(pnorm(0,mean(dv$point), mean(dv$se), lower.tail=T), pnorm(0,mean(dv$point), mean(dv$se), lower.tail=F)))

# weighted
weight_point <- sum(dv$weight*dv$point)
weight_se <- sum(dv$weight*dv$se)
weight_z <- sum(dv$weight*dv$z)

c(weight_point, weight_se, weight_z, max(pnorm(0,weight_point, weight_se, lower.tail=T), pnorm(0,weight_point, weight_se, lower.tail=F)))
Sala_i_Martin_weight[N,] <- c(weight_point, weight_se, weight_z, max(pnorm(0,weight_point, weight_se, lower.tail=T), pnorm(0,weight_point, weight_se, lower.tail=F)))




# Graph
quartz(type="pdf", width=5, height=5, file="/Users/IV/Dropbox/Nazi/Output/EBA_High_Decorated_Signal_Strength_ALL.pdf")
plot(NULL, xlim=c(1, length(dv$number)), ylim=c(min(dv$lower), max(dv$upper)), ylab=expression(paste(beta, " Radio Signal Strength")), xlab="", main="EBA - Highly Decorated")
for(i in 1:length(dv$number)){
	lines(c(i,i),c(dv$lower[i], dv$upper[i]), lwd=1.5, col="grey")
}
points(dv$number, dv$point, pch=16,cex=0.5)
abline(h=0, lwd=2, col="red")
dev.off()







## DV: Decorations: # of medals
## IV: Radio Signal Strength
N <- N+1

id <- unlist(
	lapply(1:choice.n,
		function(i)combn(1:choice.n,i,simplify=FALSE)
		)
	,recursive=FALSE)

formulas <- sapply(id, function(i)
	paste("n_medals ~ 1", sig_str, must_haves, population, region_fixed_effects, paste(choice[i],collapse=" + "), sep=" + ")
	)
	
savefunc <- function(i){
	m <- lm(as.formula(i), data= Combined_Data_Pop)
	f <- summary(m)$r.squared
	c(summary(m)$coef[2,], f)	
}

X <- lapply(formulas,savefunc)

point <- NULL
se <- NULL
z <- NULL
r2 <- NULL

for(i in 1:length(X)){
	point <- c(point, X[[i]][[1]])
	se <- c(se, X[[i]][[2]])	
	z <- c(z, X[[i]][[3]])
	r2 <- c(r2, X[[i]][[5]])		
}

lower <- point + qnorm(0.025)*se
upper <- point + qnorm(0.975)*se

dv <- data.frame(point, lower, upper, se, z, r2)
dv <- dv[order(dv$point),]
colnames(dv) <- c("point", "lower", "upper", "se", "z", "r2")
dv$number <- seq(1:dim(dv)[1])
dv$weight <- dv$r2/sum(dv$r2)

# Levine-Renelt
dv$lower[1]
dv$upper[dim(dv)[1]]

c(dv$point[1], dv$se[1], dv$z[1])
c(dv$point[dim(dv)[1]], dv$se[dim(dv)[1]], dv$z[dim(dv)[1]])

# Sala--i-Martin
# unweighted:
c(mean(dv$point), mean(dv$se), mean(dv$z), max(pnorm(0,mean(dv$point), mean(dv$se), lower.tail=T), pnorm(0,mean(dv$point), mean(dv$se), lower.tail=F)))
Sala_i_Martin_unweight[N,] <- c(mean(dv$point), mean(dv$se), mean(dv$z), max(pnorm(0,mean(dv$point), mean(dv$se), lower.tail=T), pnorm(0,mean(dv$point), mean(dv$se), lower.tail=F)))


# weighted
weight_point <- sum(dv$weight*dv$point)
weight_se <- sum(dv$weight*dv$se)
weight_z <- sum(dv$weight*dv$z)

c(weight_point, weight_se, weight_z, max(pnorm(0,weight_point, weight_se, lower.tail=T), pnorm(0,weight_point, weight_se, lower.tail=F)))
Sala_i_Martin_weight[N,] <- c(weight_point, weight_se, weight_z, max(pnorm(0,weight_point, weight_se, lower.tail=T), pnorm(0,weight_point, weight_se, lower.tail=F)))

# Graph
quartz(type="pdf", width=5, height=5, file="/Users/IV/Dropbox/Nazi/Output/EBA_Num_Decorated_Signal_Strength_ALL.pdf")
plot(NULL, xlim=c(1, length(dv$number)), ylim=c(min(dv$lower), max(dv$upper)), ylab=expression(paste(beta, " Radio Signal Strength")), xlab="", main="EBA - Number of Medals")
for(i in 1:length(dv$number)){
	lines(c(i,i),c(dv$lower[i], dv$upper[i]), lwd=1.5, col="grey")
}
points(dv$number, dv$point, pch=16,cex=0.5)
abline(h=0, lwd=2, col="red")
dev.off()



colnames(Sala_i_Martin_unweight) <- c("Beta", "SE", "Z", "CDF")
colnames(Sala_i_Martin_weight) <- c("Beta", "SE", "Z", "CDF")
Sala_i_Martin_unweight <- as.data.frame(Sala_i_Martin_unweight)
Sala_i_Martin_weight <- as.data.frame(Sala_i_Martin_weight)
Sala_i_Martin_unweight[,5] <- c("Decorations", "Highly Decorated", "# of Medals")
Sala_i_Martin_weight[,5] <- c("Decorations", "Highly Decorated", "# of Medals")


write.dta(Sala_i_Martin_unweight, "/Users/IV/Dropbox/Nazi/Output/Sala_i_Martin_unweight_ALL.dta")
write.dta(Sala_i_Martin_weight, "/Users/IV/Dropbox/Nazi/Output/Sala_i_Martin_weight_ALL.dta")










